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Abstract. A mathematical framework for optimal bilinear control of nonlin- 
ear Schrodinger equations of Gross— Pitaevskii type arising in the description 
of Bosc-Einstcin condensates is presented. The obtained results generalize 
earlier efforts found in the literature in several aspects. In particular, the 
cost induced by the physical work load over the control process is taken into 
account rather then often used L 2 - or .H^-norms for the cost of the control 
action. Well-posedness of the problem and existence of an optimal control is 
proven. In addition, the first order optimality system is rigorously derived. 
Also a numerical solution method is proposed, which is based on a Newton 
type iteration, and used to solve several coherent quantum control problems. 



1. Introduction 

1.1. Physics background. Ever since the first experimental realization of Bose- 
Einstein condensates (BECs) in 1995, the possibility to store, manipulate, and 
measure a single quantum system with extremely high precision has provided great 
stimulus in many fields of physical and mathematical research, among them quan- 
tum control theory. In the regime of dilute g BEC, consisting of N particles, 
can be modeled by the Gross-Pitaevskii equation [25] . i.e. a cubically nonlinear 
Schrodinger equation (NLS) of the form 

h 2 

ihd t i(j = --—Atp + U(x)i) + Ng\i)\ 2 ijj + W(t,xU : x € IR 3 , teR, 
2m 

with m denoting the mass of the particles, H Planck's constant, g = Anh 2 a sc /Tn, 
and a sc G K their characteristic scattering length, describing the inter-particle 
collisions. The function U (x) describes an external trapping potential which is 
necessary for the experimental realization of a BEC. Typically, U(x) is assumed 
to be a harmonic confinement. In situations where U(x) is strongly anisotropic, 
one experimentally obtains a quasi one-dimensional ( "cigar-shaped" ) , or quasi two- 
dimensional ("pancake shaped") BEC, see for instance [T7]. In the following, we 
shall assume U(x) to be fixed. The condensate is consequently manipulated via 
a time-dependent control potential W(t, x), which we shall assume to be of the 
following form: 

W(t,x) = a(t)V(x), 
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Here, a(t) denotes the control parameter (typically, a switching function acting 
within a certain time-interval [0, T]) and V(x) is a given potential. In our context, 
the potential V(x) models the spatial profile of a laser field used to manipulate the 
BEC and a(t) its intensity. 

The problem of quantum control, i.e. the coherent manipulation of quantum 
systems (in particular Bose-Einstein condensates) via external potentials W(t,x), 
has attracted considerable interest in the physics literature, cf. [H [lOj [HI [16l l24[ 
I26[ 131 j - From the mathematical point of view, quantum control problems are a 
specific example of bilinear control systems . It is known that linear or nonlinear 
Schrodinger-type equations are in general not exactly controllable in, say, L 2 (R 3 ), 
cf. [25] • Similarly, approximate controllability is known to hold for only some 
specific systems, such as [20]. More recently, however, sufficient conditions for 
approximate controllability of linear Schrodinger equations with purely discrete 
spectrum have been derived in [S]. In [T5] these conditions have been shown to be 
generically satisfied, but, to the best of our knowledge, a generalization to the case 
of nonlinear Schrodinger equations is still lacking. 

The goal of the current paper is to consider quantum control systems within 
the framework of optimal control, cf. |29j for a general introduction, from a partial 
differential equation constrained point of view. The objective of the control pro- 
cess is thereby quantified through an objective functional J = J(tp,a), which is 
minimized subject to the condition that the time-evolution of the quantum state 
is governed by the Gross-Pitaevskii equation. Such objective functionals J{ij), a) 
usually consist of two parts, one being the desired physical quantity (observable) 
to be minimized, the other one describing the cost it takes to obtain the desired 
outcome through the control process. In quantum mechanics, the wave function 
if>(t, ■) itself is not a physical observable. Rather, one considers self-adjoint linear 
operators A acting on i/)(i, •) and aims for a prescribed expectation value of A at 
time t — T > 0, the final time of the control process. Such expectation values are 
computed by taking the L 2 -inner product (ip(T, -),Aip(T, -)).L 2 (R d )- Note that this 
implies that the corresponding ip(t, •) is only determined up to a constant phase. 
This fact makes quantum control less "rigid" when compared to classical control 
problems in which one usually aims to optimize for a prescribed target state. 

There are many possible ways of modeling the cost it takes to reach a certain 
prescribed expectation value. The corresponding cost terms within J(ip, a) are 
often given by the norm of the control a(t) in some function space. Typical choices 
are L 2 (Q,T) or iJ 1 (0,T). However, these choices of function spaces for a(t) often 
lack a clear physical interpretation. In addition, cost terms based on, say, the L 2 - 
norm of a tend to yield highly oscillatory optimal controls due to the oscillatory 
nature of the underlying (nonlinear) Schrodinger equation. The same is true for 
quantum control via so-called Lyapunov tracking methods, see, e.g., |12j . In the 
present work we shall present a novel choice for the cost term, which is based on 
the corresponding physical work performed throughout the control process. 

We continue this introductory section by describing the mathematical setting in 
more detail. 

1.2. Mathematical setting. We consider a quantum mechanical system described 
by a wave function ip(t, ■) <G L 2 (R d ) within d = 1, 2,3 spatial dimensions. The case 
d = 1,2 models the effective dynamics within strongly anisotropic potentials (re- 
sulting in a quasi one or two-dimensional BEC). The time-evolution of ip(t,-) is 
governed by the following generalized Gross-Pitaevskii equation (rescaled into di- 
mensionless form): 

(1.1) idt^ = -^Atlj + U(x)iP + X\^\ 2a iP + a(t)V(x)i;, xeR d ,teR, 
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with A ^ 0, (J < 2/(d — 2), and subject to initial data 

ip{0, •) = ip £ L 2 (R d ), a(0) = a G E. 
For physical reasons we normalize HV'ollz, 2 ^) = lj which is henceforth preserved 



by the time-evolution of (1.1). In addition, the control potential is assumed to be 



V G W 1 '°°(R d ), whereas for U(x) we require 

U G a°°(m d ) such that d k U £ L°°{R d ) for all multi-indices k with \k\ > 2. 
In other words, the external potential is assumed to be smooth and subquadratic. 



One of the most important examples is the harmonic oscillator U(x) — ^\x\ 2 . Due 
to the presence of a subquadratic potential, we restrict ourselves to initial data ipo 
in the energy space 

(1.2) E :={?/> € H 1 (R d ) : xip £ L 2 (R d )} . 

In particular, this definition guarantees that the quantum mechanical energy func- 
tional 

(1.3) E(t) = f hv^{t 1 x)\ 2 + ^-m i x)\^ +2 + {a{t)V(x) + U{x))m,x)\ 2 dx 

J R d I a + 1 

associated to ( |1.1[ ) is well defined. 

Remark 1.1. Note that a < 2/{d — 2) allows for general power law nonlinearities 
in dimensions d — 1,2, whereas in d = 3 the nonlinearity is assumed to be less 
than quintic. From the physics point of view a cubic nonlinearity a = 1 is the 
most natural choice, but higher order nonlinearities also arise in systems with more 
complicated inter-particle interactions, in particular in lower dimensions; compare 
|17j . From the mathematical point of view, it is well known that the restriction 
(7 < 2/(d—2) guarantees well-posedness of the initial value problem in the energy 
space E; see [7[ 0. In addition, the condition A ^ (defocusing nonlinearity) 



guarantees the existence of global in-time solutions to (1.1 1; see [7]. Hence, we do 
not encounter the problem of finite-time blow-up in our work. 

Although ( |1.1[ ) conserves mass, i.e. \\tp(t, •)||z, 2 (R<i) = I! V'oIIl 2 ^) f° r an £ G R, the 
energy E(t) is not conserved. This is in contrast to the case of time-independent 
potentials. In our case, rather one finds that 

(1.4) ~E(t)=a(t) f V(x)m,x)\ 2 dx. 
at j R d 

The physical work performed by the system within a given time-interval [0, T] is 
therefore equal to 

(1.5) E{T) - E(0) = [ a(t) [ V(x)\ip(t,x)\ 2 dx dt. 

JO JR d 

Thus a control a(t) acting for t £ [0, T] upon a system described by ( 1.1 ) requires a 



certain amount of energy, which is given by ( 1.5 ). It, thus, seems natural to include 
such a term in the cost functional of our problem in order to quantify the control 
action. 

Indeed, for any given final control time T > 0, and parameters 71 ^ 0, 72 > 0, 
we define the following objective functional: 

(1.6) J(V,a) := (i>(T,-),A*P(T,-))l 2m + 7l ^ (E(t)) 2 dt + 72 J (a(t)) 2 dt, 

where A : E — > L 2 (R d ) is a bounded linear operator which is assumed to be essen- 
tially self-adjoint on L 2 (R d ). In other words, A represents a physical observable 
with spec (A) C R. A typical choice for A would be A = A' — a where a £ R is 
some prescribed expectation value for the observable A' in the state ip(T,x). For 
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example, if a G spec (A') is chosen to be an eigenvalue of A' , the first term in 
J(ijj, a) is zero as soon as the target state ip(T, •) is, up to a phase factor, given by 
an associated eigenfunction of A' . However, one may consider choosing a G K such 
that it "forces" the functional to equidistribute between, say, two eigenfunctions. 

Remark 1.2. We also remark that in the case A = P„ — 1, where P v denotes the 
orthogonal projection onto a given target state ip e L 2 (R d ), the first term on the 



right hand side of (1.6) reads 



(1.7) WT,-),A^(T,-)) LHRd) = K^.^OW*)! 2 -!, 



using the fact that \\ip(T, •) \\ L 2 / ffi <j) = 1. Expression (1.7) is the same as used in 
recent works in the physics literature; see |14j . 



Using (1.4), we find that the objective functional J(ip,a) explicitly reads 



J(i>,a) := (iP(T,-),AiP(T,-))l Hm 

+ 7l f {a{t) f ( f V(x)\i>(t,x)\ 2 dx) dt + 7 2 f (d(t)) a dt. 

JO \JR d J JO 

Here, the second line on the right hand side displays two cost (or penalization) 
terms for the control: The first one, involving 71 0, is given by the square of 



the physical work, i.e. the right hand side of (1.4). The second is a classical cost 
term as used in [14] . In our case, the second term is required as a mathematical 
regularization of the optimal control problem, since for general (sign changing) 
potentials V G L°°(R d ) the weight factor 

(1.9) u(t) := [ V{x)\iP(t,x)\ 2 dx 

might vanish for some t£l. In such a situation, the boundedness of variations of 
a(t) is in jeopardy and the optimal control problem lacks well-posedness. Hence, 
we require 72 > for our mathematical analysis, but typically take 72 <C 71 in our 
numerics in Section [5] to keep its influence small. Note, however, that in the case 
where the control potential satisfies the positivity condition 

V{x) > 5 > Vx G R d , 

we may choose 72 = and all of our results remain valid. 

Remark 1.3. In situations where the above positivity condition on V(x) does not 
hold, one might think of performing a time-dependent gauge transform of ip, i.e. 

^{t : x^j = exp f — xk J a(s) ds^Jip(t, x), 

with a constant k > min x ^dV(x), assuming that the minimum exists. This yields 
a Gross-Pitaevskii equation for the wave function ip with modified control potential 
V(x) = (k + V(x)) > for all x G M d . Note, however, that this gauge transform 



leaves the expression (1.8) unchanged and hence does not improve the stuation. 
Only if one also changes the potential V(x) within J(ip,a) into V(x), the problem 
does not require any regularization term (proportional to 72). Note, however, that 
such a modification yields a control system which is no longer (mathematically) 
equivalent to the original problem. In fact, replacing V{x) by V(x) in the objective 
functional J(if),a) corresponds to increasing the parameter 72 by k. 
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1.3. Relation to other works and organization of the paper. The math- 
ematical research field of optimal bilinear control of systems governed by partial 
differential equations is by now classical, cf. [131 118) for a general overview. Surpris- 
ingly, rigorous mathematical work on optimal (bilinear) control of quantum systems 
appears very limited, despite the physical significance of the involved applications 
(cf. the references given above). Results on simplified situations, as, e.g., for finite 
dimensional quantum systems, can be found in [5] (see also the references therein). 
More recently, optimal control problems for linear Schrodinger equations have been 
studied in [21 [T5], In addition, numerical questions related to quantum control 
are studied in 3 , 30] . Among these papers, the work in [15] appears closest to our 
effort. Indeed, in |15j . the authors provide a framework for bilinear optimal control 
of abstract (linear) Schrodinger equations. The considered objective functional in- 
volves a cost term proportional to the L 2 -norm of the control parameter a(t). The 
present work goes beyond the results obtained in 15J in several repects: First, we 
generalize the cost functional to account for oscillations in a(t) and in particular 
for the physical work load performed throughout the control process. In addition, 
we allow for observables A which are unbounded operators on L 2 . Second, we 
consider nonlinear Schrodinger equations of Gross-Pitaevskii type, including un- 
bounded (subquadratic) potentials, which are highly significant in the quantum 
control of BECs. This type of equation makes the study of the associated control 
problem considerably more involved from a mathematical point of view. 

The rest of this work is organized as follows. In section [2] we clarify existence of a 
minimizer for our control problem. In particular, we prove that the corresponding 
optimal solution tfj^(t,x) is indeed a mild (and not only a weak) solution of ( |1.1[ ), 
depending continuously on the initial data ipo- Then, in section [3] the adjoint 
equation is derived and analyzed with respect to existence and uniqueness of a 
solution. It is our primary tool for the description of the derivative of the objective 
function reduced onto the control space through considering the solution of the 
Gross-Pitaevskii equation as a function of the control variable a. The results of 
section [3] are paramount for the derivation of the first order optimality system in 
section [4| In section [5] a gradient- and a Newton- type descent method are defined, 
respectively, and then used for computing numerical solutions for several illustrative 
quantum control problems. In particular, we consider the optimal shifting of a linear 
wave package, splitting of a linear wave package and splitting of a BEC. The paper 
ends with conclusions on our findings in section [6] 

Notation. Throughout this work we shall denote strong convergence of a se- 
quence (x n ) I j S N by x n — > x and weak convergence by x n — 1 x. For simplicity, we 
shall often write ip(t) = ip(t, •) and also use the shorthand notation L\ L q x instead of 
L p (0,T;L«(M d )). Similarly, H\ stands for H 1 ^^), with dual Hf 1 = (F x (0,T))*. 

2. Existence of minimizers 

We start by specifying the basic functional analytic framework. For any given 
T > 0, we consider if 1 (0, T) as the real vector space of control parameters a(t) G M. 
It is known [5] that for every a G H (0,T), there exists a unique mild solution 
ijj G C([0,T]; S) of the Gross-Pitaevskii equation. More precisely, ip solves 

rP(t, x) = S(t)Mx) -if S(t-s) (A|V(s, -)\ 2a ^(s, •) + a{s)Vij(s, •)) (x)ds, 
Jo 

where from now on we denote by 



(2.1) 



S(t) = e-* tH , H = --A + U(x), 
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the group of unitary operators {S(t)}t£R generated by the Hamiltonian H. In other 
words, S(t) describes the time-evolution of the linear, uncontrolled system. Next, 
we define 

(2.2) T(0,T) :=L a (0,T;E)n.ff 1 (0,T;E*), 

where E* is the dual of the energy space S. Then the appropriate space for our 
minimization problem is 

A(0,T) := {(ip,a) G T(0,T) x H\0,T) : tp is a mild solution of ((TlJ }. 

Since the control a is real-valued, it is natural to consider A(0,T) as a real vector 
space and we shall henceforth equip L 2 (R d ) with the scalar product 



(2.3) (£»L2(Rd) = Re / £(x)ip(x) dx, 

which is subsequently inherited by all L 2 -based Sobolev spaces. (Note that this 
choice is also used in [8].) From what is said above, we infer that the space A(0, T) 
is indeed nonempty. 

With these definitions at hand, the optimal control problem under investigation 
is to find 

(2.4) J,= inf J(^,a). 

(V,a)eA(0,T) 

We are now in the position to state the first main result of this work. 

Theorem 2.1. Let A > 0, < a < 2/{d-2), V G W^°°(R d ), and U G C°°{R d ) be 
subquadratic. Then, for any T > 0, any initial data ?/>o G S, «o G K and any choice 



of parameters 71 > 0, 72 > the optimal control problem (2.4) has a minimizer 
(i,a,)€A(0,T). 



The proof of this theorem will be split into three steps: In subsection |2.1| we 
shall first prove a convergence result for minimizing, or more precisely, infrmizing 



sequences. We consequently deduce in subsection 2.2 that the obtained limit -0* 



is indeed a mild solution of (1.1). Finally, we shall prove lower semicontinuity of 



J(4>, a) with respect to the convergence obtained before. 

2.1. Convergence of infimizing sequences. First note that there exists at least 
one infimizing sequence with an infimum —00 J* < +00, since A(0,T) and 
J : A(0, T) — » R. Then we have the following result for any infimizing sequence. 

Proposition 2.2. Let (ip n , a n )neN 06 an infimizing sequence of the optimal control 



problem given by (1.6). Then under the assumptions of Theorem \ 2. 1\ there exist 
a subsequence, still denoted by (^ n ,a n ) n gn, and functions a* G H L (0, T), ip* G 
L°°(0,T;S), such that 

a n — a* in H 1 (0, T) , and a n —> a* in L 2 (0, T) , 

ipn "0* in £ 2 (0,T; S), 

V>„ -> ^* in L 2 (0, T; L 2 (R d )) n L 2 (0, T; L 2a+2 (R d )), 
as n — > +00. Furthermore it holds that 

(2.5) ip n (t) -> in i 2 (M d ), and ^ n (i) ip*(t) m E 

/or almost all t G [0, T]. 

Proof. By definition, J and thus it is bounded from below. For an infimizing 
sequence (^„,a„)„gM the sequence of objective functional values ( J(ip n , a n ))neN 
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converges and is bounded on K. Hence, it holds that J(i/j n , <*n) ^ C < +00 for all 
n G N. Since 72 > it follows that 

i-T 

(a n (t)) 2 dt^C < +00. 

For smooth a n : [0, T] — > R we compute 

t / t \ ^ 2 

a„(i) = a„(0) + y d„(s) (is a„(0) + i T J (««( s )) 2 rfs j < +°°' 

and thus a„ is bounded in L°° (0, T) . By approximation (using the fact that a n (0) = 
«o is fixed), the sequence (a n ) n ^ is uniformly bounded in L°°(0, T), which in turn 
implies a uniform bound in L 2 (Q,T) and thus in i? 1 (0,T). Hence, there exists a 
subsequence, still denoted (a„) nS N, and a* € i? 1 (0,T), such that 

a„ a, G ff^OjT). 

Moreover, since i? 1 (0, T) is compactly embedded into L 2 (0,T), we deduce that 
a n — > a* in L 2 (0,T). Next, we recall that 

d f 

j t E n {t) = a n (t) J ^ V(x)\i) n (t,x)\ 2 dx 

and hence 

\\E n \\ L % < KIU?IMU?Wl£s' 
in view of mass conservation ||'0„(<)|| i 2 = ||V>o||z, 2 • Since E n (0) = E depends only 
on %pQ and ao (and is thus independent of n £ N), the same argument as before 
yields ||^ n ||i« ^ C. Recalling the definition of the energy (1.3 1 and the fact that 
A ^ 0, we obtain 

(2.6) ^||W>„(i)|| 2 2 ^ \\E n \\ Lr +c||a n || i oc||^ ||| 2 + C\\xip n (t)\\h 
z 



again using conservation of mass HV'nWIli 2 = llV'olli 2 - Furthermore, it holds that 
' \x\ 2 \ip\ 2 dx = 2Re / i\x\ 2 ^ ( ^Aip - X\tp\ 2,J ip - XaVip - Uifi j dx 



dt .lOd .lOd \ 2 

^ M_././j.Ml2 , llw././Vil, 

|Z 2 > 



21m / xi^Vip dx < \\xip{t)\\ 2 L2 + ||V^(£)|| 



which, in view of the bound (2.6 1 and Gronwall's inequality yields 



^c(\\E n \\ Lr + \\a n \\ Lr \\M\] 
for all t £ [0, T]. In summary we have shown 

(2-7) = \\Mt)\\ 2 Hi + ||a^n(*)lli; < C 

where C > is independent of n £ N and t £ [0, T}. Hence, ip n is uniformly bounded 
in i°°(0,T;S) and in particular in L 2 (0,T;£). By reflexivity of L 2 (0,T;£), we 
consequently infer the existence of a subsequence (denoted by the same symbol) 
such that 

ij) n — -0* in L 2 (0, T; £) as n —> +00. 

To obtain the strong convergence announced above, we first note that ( |1.1| implies 
d t ipn G L°°(0,T;£*). On the other hand, £ is compactly embedded in L 2 (R d ). 
Thus, we can apply the Aubin-Lions Lemma to deduce 

ipn V>* in L 2 ((0,T) x R d ). 
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In particular, there exists yet another subsequence (still denoted by the same 
symbol), such that 

ipn{t) ^^-^ m L 2 (R d ), for almost all t € [0,T]. 

In order to obtain weak convergence in the energy space, i.e. ip n (t) ~^ i^*{t) in E, 
we fix t £ [0,T] such that ip n (t) -> ip*(t) in L 2 (R d ). In view of 07f\ , every 
subsequence of if> n (t) has yet another subsequence such that ip n (t) converges weakly 
in E to some limit. On the other hand, this limit is necessarily given by ip*(t), since 
ip n {t) -> m L 2 (R d ). Hence the whole sequence converges weakly in E to 

ij)*(t). By lower-semicontinuity of the E-norm we can deduce ||?/Vi(i)||£ S% C and 
thus ip* G L°°(0,T;E). 

Finally, the announced convergence in L 2 (0, T; L 2a+2 (R d )) is obtained by invok- 
ing the Gagliardo-Nirenberg inequality, i.e. 

(2-8) llellL^C||e||^ M ||Ve||i ( 2 r) , 



where 2 < r < and 5(r) — d{\ — ^). This concludes the proof of Proposition 



2.2 



2.2. Minimizers as mild solutions. Next we prove that the limit ip* obtained in 



the previous subsection is indeed a mild solution of ( 1.1 ) with corresponding control 
a*. From the physical point of view, this is important since it implies continuous (in 
time) dependence of ip* upon a given initial data ipQ. To this end, one should also 
note that H 1 (0,T) <^-» C(0,T) (using Sobolev imbcddings), and hence the obtained 
optimal control parameter a*(t) is indeed a continuous function on [0,T]. 

Proposition 2.3. Let (-0*, a*) £ T(0, T) x H 1 (0, T) 6e t/ie limit obtained in Propo- 



sition 2.2 Then is a mild solution of (1.1 1 with control a* and 



^ £ C([0,T];E)nC 1 ([0,T];E*). 



In particular, this implies that the convergence result (2.5) holds for all t £ [0, T]. 
Proof. First we note that, by construction, each ip n satisfies 

ij} n {t) = S{t)i> -i f S(t-s)(X\Ms)\ 2a Ms)+a n (s)VMs))ds 



t 

2<r„ 

oyu — s) y^yifJnya j\ 



for all t £ [0, T]. Here and in the following we shall suppress the x-dependence of 
ijj for notational convenience. In order to prove that ip* is a mild solution corre- 
sponding to the control a* , we take the L 2 -scalar product of the above equation 
with a test function x G C£°(]R d ). This yields 

(ipn(t),x)Ll = (S{t)Tjj , X ) L 2-i\ [ (S{t-s)\^ n (s)\ 2 ^ n {s), X ) L *ds 

(2.9) Jo 

-i / (S(t - s)a n (s)Vip n (s),x) L 2 ds. 
Jo 

In view of Proposition 2.2 the term on the left hand side of this identity converges 
to the desired expression for almost all t £ [0,T], i.e. 

lim (ip n (t),x)Ll = (i>*(t),x)L*- 

n— >oo x x 

In order to proceed further, we note that for any / 6 T>'(W l ) it holds that 

(2.10) (S(t-s)f(s),xhl = (f(s),S(s-t) X )Ll, 
and we therefore define 

(2.11) rMlxK^C, x^x(;x)~S(--t) X (x), 
for which we can prove the following regularity properties. 
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Lemma 2.4. There exists a constant C = C(T) > such that for all t € [0, T] it 
holds that 

sup (\\xx(s)\\ L 2 + ||Vx(s)|U|) < C(T) < +oo, 
se[o,t] 



where the function x is defined in (2.11). In particular, the function x is bounded 
in L°°(0,t;L 2<T +' 2 (R d )). 



Proof of Lemma \2.4\ The norm ||x( s )l|i 2 = \\S( S ~ 1 )x\\l 2 is conserved since S(t) 
is a unitary operator on L 2 (R d ). Furthermore, it holds that 

idt\y, S(t)} = H[W, S(t)] + [V, H]S(t) = H[W, S(t)] + Vf/S(t), 

and hence 

[V, S(t)} = -i [ S(t - s)VU S(s) ds. 
Jo 



We can thus estimate 

l|Vx(*)lk - \\VS(t-s)x\\ L * 



^ \\S{t - s)v x \\li + 

r t-s 



S(t-s- t)VUx(t) dr 



(2.12) 

^\\Vxhl+C I \\xx(t)\\ L 2 dr, 
Jo 

since U is subquadratic, i.e. |VtT(x)| ^ C\x\. Likewise, we deduce 

[x, S(t)] = -i [ S(t - s)VS(s) ds 
Jo 

and hence 



LI 



(2.13) 



\\ x x{s)\\li < Wxxhi 



\^x(t)\\li dr. 



Combining the estimates ( 2.12 ) and ( 2.13 1 and applying Gronwall's inequality yields 

\\xx{s)\\li + \\Vx(s)hl ^C(\\x X \\li + ||Vx|U 2 ) < +oo, 

where C > 0. The bound in L x (0,t; L 2<T+2 (R rf )) then follows fr om the uniform-in- 
time bound in H 1 (M. d ) and the Gagliardo-Nirenberg inequality (2.8). □ 



With the result of Lemma 2.4 at hand, we consider the second term on the right 



hand side of ( |2.9[ ). Rewriting it using (2.10), we estimate 
' (\Ms)\ 2a Ms)-\^(s)\ 2a Ms),x(s)) Ll ds 

\ip n (s, x)\ 2a ipn(s, x) - \i/j sr (s,x)\ 2a ip*(s,x)\ \x(s,x)\ dx ds 

(\ip n (s,x)\ 2cr + \ip*(s 1 x)\ 2a ) \ip n (s,x) - ^(s,x)\\x(s,x)\ dxds. 



< C 



lo Jm. d 

By Holder's inequality, it holds that 

|2<7 , /„ „m2ct 



JR d 



\ip n (s,x)\ a + |V>*(s,a;)| a \\ipn(s, x) - ip*(s,x)\\x(s)\ dxds 



\ t x t x / t t 



where, in view of Lemma 



2.4 



we have ||x|| i oo i 2<7+2 < +oo. In addition, Proposi- 
tion [2?2] implies that the factor inside the parentheses is bounded and that 

lim \\ip n — ^ >t || L 2 L 2^+2 =0. 
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Thus, we have shown that the second term on the right hand side of (2.9) vanishes 
in the limit n ~ > oo. 

It remains to treat the last term on the right hand side of (2.9), rewritten via 



(2.10). We first estimate 

rt 

(a n (s)Vip n (s) - a*(s)Vip*(s),x(s)) L2 ds 



^11 \a n (s)\\V(x)\\ip n (s,x)-il}*(s,x)\\x(s,x)\dxds 
\a n (s) - a*(s)| |V(a;)| |^*(s,a;)||x(s,a;)| dxds. 



Here, the last term on the right hand side can be bounded by 
r t t 

\a n (s) - a*(a)\\V(x)\\ip*(8,x)\\x(s,x)\ dxds 

«S \\a n - a*|U=||V||i«»||V'|U|i|||xlUr'il °j 

in view of the convergence of a n — > a* in L 2 (0,T). For the remaining term we use 
the fact that V £ L°°(M. d ) and Holder's inequality to obtain that 
r t t 

\a n (s)\ \V(x)\ \ip n (s,x) - ip*(s,x)\\x(s,x)\ dxds 

^ ll"n|Uf||^|Us°||V'n -^*llifilllxlUfi» "^°°> P. 

due to the results of Proposition |2.2| and Lemma |2.4| 

In summary this proves that tp* G T(0,T) satisfies, for almost all t £ [0, T], 

M*) = S(t)^>o -if S{t~s) (\\^(s)\ 2a ^(s) + a«(«)VV.(*)) ds, 
Jo 

i.e. "0* is a weak T,-solution in the terminology of [SI Definition 3.1.1] (where the 
analogous notion of weak i? 1 -solutions is introduced). In order to obtain that ip* 
is indeed a mild solution we note that 

V>* € T(0,T) C{[0,T}; L 2 (M. d )) D C{[0,T}; L 2a+2 {R d )) 



by interpolation and the Gagliardo-Nirenberg inequality ( |2.8[ ). Classical arguments 
based on Strichartz estimates then yield uniqueness of the weak S-solution ip*. 
Arguing as in the proof of [8, Theorem 3.3.9], we infer that ip* is indeed a mild 
solution to (jOJ, satisfying 0* € C([0,T];S)nC rl ([0,T];E*). □ 

2.3. Lower semicontinuity of objective functional. In order to conclude that 
the pair (if}*, a*) G A(0, T) is indeed a minimizer of our optimal control problem, it 
remains to show lower semicontinuity of the functional J(ip, a) with respect to the 
convergence results established in Proposition |2.2| 



Lemma 2.5. For the sequence constructed in Proposition \2.£\ it holds that 
J* = liminf J(ip n ,a n ) > J(ip*, a*)- 

n— >oc 

Proof. Since A 6 £(£, L 2 (R d )) by assumption, the sequence (A^j n (T)) n ^ con- 
verges weakly to Aip(T) in L 2 (R d ). In addition ip n (T) -> ip*(T) in L 2 (R d ) as 
7i — y oo by Proposition |2.3| and hence the estimate 



| (*Pn(T) , Aip n (T)) L 2 - (MT),AMT))li\ 

s? I (T) - ^ (T) , A^ n (T)) Li \ + \ (T) , A^ n (T) - ^* (T) )) l; I 
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yields convergence of the corresponding term in the objective functional ( 1.8 1. Next, 
we consider the cost term involving 71. In view of (1.9), we define 

V(x)\Mt,x)\ 2 dx, w«(t):= 



and estimate 



V{x)\^{t,x)\ 2 dx, 



lim inf 

n— voo 



(d„(i))X(t)d* ^ 



liminf / (d n (t)) 2 uj 2 (t) dt + liminf / {d n {t)f {u 2 n {t) - u; 2 (t)) dt 



(2.14) 



Note that < u n (t) ^ \\ V\\ L ^ \\ip \\ 2 L 2 independently of n e N and t e [0,T] and 
that the same holds for uj*(t). The first term on the right hand side of (2.14) is 
convex in a„ and thus satisfies 



(2.15) 



lim inf 

n— 



(a n {t)) 2 tu 2 {t) dt ^ / (d«{t)) 2 uj 2 (t) dt 



since any convex and lower semicontinuous functional is weakly lower semicontin- 
uous. On the other hand, Proposition |2.2| implies 

(2.16) 



lim inf oj n (t) ^ oj* (t) ^ for all t € [0, T]. 



Thus, using (2.15) and (2.16) together with Fatou's Lemma yields 



lim inf 

n— f 00 



(a n (t))W n (t) dt 



(d4t)) 2 cu 2 (t) dt + / lim inf (d„(i)) 2 liminf (wl{t) - ul{t)) dt 



(d*(i)) 2 w 2 (i) dt. 



Finally the cost term involving 72 is lower semicontinuous by convexity and weak 
convergence of a n in if 1 (0,T). □ 

In summary, we have shown that J* = lim inf ^^oo J(ip n , a n ) ^ J(V>*,a*) and 
thus indeed J* = J(?A* , a*). In other words, (")/>*, a*) € A(0,T) solves the optimiza- 
tion problem. 

Remark 2.6. Note that the bound on xijj n (t, •) in L 2 (M. d ), obtained in Proposition 



2.2 is indeed crucial for proving the weak lower-semicontinuity of J(i/j, a). Without 



such a bound on the second moment, we would only have 

Mt)^^m ™L 2 oc (R d ), 

due to the lack of compactness of H 1 (R d ) L 2 (M. d ). In this case, the lower 
semi-continuity of the term (ip(T), Aip(T)) L 2 is not guaranteed. A possible way 
to circumvent this problem would be to assume that A is positive definite, which, 
however, is not true for general observables of the form A = A' — a, with a € M. A 
second possibility would be to assume that A is localizing, i.e. for all i/j £ H 1 (M. d ): 
s\rpTp xem d(Aip(x)) C B(R), for some R < +00. 
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3. Derivation and analysis of the adjoint equation 

In order to give a characterization of a minimizer (-0*, a*) € A(0, T), we need to 
derive the first order optimality conditions for our optimal control problem (2.4). 
For this purpose, we shall first formally compute the derivative of the objective 
functional J(ip, a) in the next subsection and consequently analyze the resulting 
adjoint problem. A rigorous justification for the derivative will be given in Section 

SI 

3.1. Identification of the derivative of J(ip, a). The mild solution of the nonlin- 
ear Schrodinger equation (1.1), corresponding to the control a € iP(0,T), induces 
a map 

ip : iI 1 (0,T) ->■ T(0,T) : a i-» 0(a). 

Using this map we introduce the unconstrained or reduced functional 

J : iP^(0,T) ->R, a^- J(a) := J{ij){a),a). 

For the characterization of critical points, we need to compute the derivative of 
J. For this calculation let 5 a £ i? 1 (0,T) with S a (0) — be a feasible control 
perturbation. (Recall that -ff 1 (0, T) <-» C(0, T) and hence it makes sense to evaluate 
S a (t) at t = 0.) Then the chain rule yields 

(j'(a),S a ) H -i H i = (d^,J(ip(a),a),ip'(a)S a ) r * r 

(3.1) 

where T* denotes the dual space of T = T(0, T) for any given T > 0. The main dif- 
ficulty lies in computing ip'( a ) since is given only implicitly through the nonlinear 
Schrodinger equation (1.1). 

In the following, we shall write the (nonlinear) partial differential equation (1.1) 
in a more abstract form, i.e. 

(3.2) P(tp, a) := id t ip -Hi/}- a{t)V{x)^ - A|0| 2<T i/; = 0, 

where H = — |A + U(x) denotes the linear, uncontrolled Hamiltonian operator. 
Setting if> — ip(a) and differentiating with respect to a formally yields 

4-P(il>(a),a) =d^P{i){a),a)ip'{a) + d a P{tl}{a),a) = 0. 
da 

Next, assuming that d^P is invertible, we solve for ip'(a) via 
if/ (a) = —d^P(ip, a)~ 1 d a P('ip(a),a). 

Thus it holds that 

(d^J(ip(a),a),i>'(a)S a ) r , T 

= {—d^J '(ip(a) , a), d^P(i>(a) , a) -1 d a P(ip(a) , a)S a ) T , r , 

which can be rewritten as 

(d-dj J(tp(a), a), ip'(a)S a ) r , T 

(3.3) 

= {-d a P{i>{a),a)*d^P{^{a),a)~*d^J{^{a),a),8 a ) H - l Hl . 
Here we abbreviate 

a^P(V»(a),a)- := (9^P(^(a), a)*)" 1 = (tyP^a), a)" 1 )*. 



Substituting (3.3) into equation ( |3.l[ ), we see that critical points of (2.4) satisfy 

= (J'(a),5 a ) H -i „i = (d a J(ip(a),a),5 a ) H -i H i + 
(3.4) ' ' 4 _ * ' * 

(—d a P(ip(a), a)*d t pP(tp(a) 7 a) *d^J(tp(a),a), S a ) H -i H i 
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for all S a £ i? 1 (0, T) such that 5 a (0) = 0. In order to obtain (3.4 1 in a more explicit 
form, we (formally) compute the derivative 

(3.5) 0*Pty, a)£ = idtt -H£- a(]t)V{x)t - A(a + 1)|V| 2 ^ - Ac#r~V?, 
acting on £ <E L 2 (M. d ) C £*. Analogously, we find 

d a P{ip,a) = -V(x)ip. 

Next, we define 

(3.6) y> := d^P(tp(a),a)~*d^J(rp(a),a) ) 



which, in view of (3.4 1, allows us to express J'ipt) £ (i/ 1 (0,T))* in the following 
form: 

(3.7) J' (a) = d a J{^{a), a) - d a P(t[>(a),a)*tp. 

We consequently obtain J 1 {a) by explicitly calculating the right hand side of this 
equation (given in (4.3) below), provided we can determine tp. 

In order to perform this calculation, we recall that the duality pairing between 
£ £ L 2 (R d ) C S* and tp £ £ can be expressed by the inner product defined in (2.3 1. 
Thus, (3.6) implies 

(3.8) (ip,d i ,P(ip(a),a)5 4 ,) L 2 L 2 = (d^J(ip(a), a), S^) L 2 L 2, 

for all test functions 5^ £ T(0,T) such that 5^,(0) = 0. This is the correct "tangent 
space" for ip in view of the Cauchy data 

4>(0) + fy(0) = i/;o and ^(0) = </V 



By virtue of the symmetry of the linearized operator d^,P{^){a), a), equation (3.8) 
corresponds to the weak formulation of the following adjoint equation: 



2- _ x_l„/.|2 ff -2„,2^ _ SJ(lp,a) 



(3.9) 



id t tp -Hip- a(t)V(x)tp - \{a + if - MM 

for all t £ [0, T] and with data: tp(T) = i 



SJ(ip, a) 



Here, SJ g^f^ denotes the first variation of J(tp, a) with respect to the value of 
tp(t) £ H 1 (M. d ), where -0 is the solution of ( |1.1[ ) with control a. Likewise, ^^fey 
denotes the first variation with respect to solutions of (1.1) evaluated at the final 
time t — T. Explicitly, these derivatives are given by 

5J{il>,a) = 4 ( d ^ 2 f I" v(x)\ip(t,x)\ 2 dx] V{x)i){t,x) 



(3.10) H{t) 

= 4(a(i)) 2 w(t)\/(a;)0(t,a;), 



in view of the definition ( 1.9 1, and 

(3.11) S -^^=mT,-),A<l>(T,.)) Li Ai>(T,x). 



The system (3.9) consequently defines a Cauchy problem for tp with data given at 
t = T, the final time. Thus, one needs to solve (3.9) backwards in time, a common 
feature of adjoint systems for time-dependent phenomena. 

Remark 3.1. In fact, tp can also be seen as a Lagrange multiplier within the 
Lagrangian formulation of the optimal control problem. In oder to see this, one 
defines the Lagrangian 



L(ip,a,tp) = J(ip,a) - (tp,P(ip,a)) 



14 M. HINTERMULLER, D. MARAHRENS, P. A. MARKOWICH, AND C. SPARBER 



where P(ij),ot) is the nonlinear Schrodinger equation given in (3.2). Formally, the 
Euler-Lagrange equations associated to L(tfj, a,p) yield (3.7) and (3.9). In Section[5] 
we shall use the Lagrangian formulation to formally compute the Hessian of the 
reduced objective functional J {a). 



In the next subsection, we shall set up an existence theory for (3.9 
turn will be used to rigorously justify the above derivation in Section |4 



, which in 
below. 



3.2. Existence of solutions to the adjoint equation. In order to obtain exis- 



tence of solutions to (|3.9|), we need sufficiently high regularity of if, the solution of 

we define 



for all multi-indices j and k with 



the Gross-Pitaevskii equation (1.1). For this purpose, for every m G 
E m := {ip € L 2 {R d ) : x j V k ip € L 2 (I 
|i| + |fe| s^m}, 
equipped with the norm (note that E 1 = E): 

Helium := E Ik'vVl 



|i|+|fc|<m 



Remark 3.2. If the external potential U(x) were in L°°{R d ), it would be enough 
to work in the space H m (R d ) instead of E m . In the presence of an external sub- 
quadratic potential, however, we also require control of higher moments of the wave 
function if with respect to x. 



Lemma 3.3. Let S(t) be given by f2jl with U G C°°(Ifr;R) and subquadratic. 



Then, there exists a constant c > such that 

\\S{t)i> \\^^e ct U4^, 
for all t G [0, oo) and ip G E m . 

The proof of this lemma can be deduced by differentiating the E m -norm with 
respect to time and applying Gronwall's inequality. It consequently implies the 



following regularity result for solutions to (1.1) 



Lemma 3.4. Let A ^ 0, a G N with a < 2/(d - 2), and U G C°°{R d ) be sub- 
quadratic. For m > d/2, let i[) Q G E m , and V G W m >°°{R d ). Then the mild 



solution of (1.1) satisfies ip G £°°(0, T; E m ). 



In view of Lemma |3.3[ this result can be proved by following the same argu- 
ments as in [SI Theorem 5.5.1]. Having obtained ip G L°°(0,T;E m ), we infer 
if) G L°°((0,T) x R d ) by the Sobolev embedding H m {R d ) ^ L°°(R d ) whenever 
m > d/2. Thus, all the '(/'"dependent coefficients appearing in adjoint equation 



(3.9) are indeed in L° 



Remark 3.5. Note that Lemma [3.4| requires us to impose a G N, which together 
with the condition a < 2/(d — 2) necessarily implies d ^ 3. The reason is that 
for general a > (not necessarily an integer) the nonlinearity \ij)\ 2 ' T, 4' is not locally 
Lipschitz in E m (cf. Lemma 4.2) and the life-span of solution ip(t,-) G E TO is in 



general not known, see [5j for more details. 

From now on, we shall always assume that V G W m,oa (R d ) for m > d/2 and 
U G C°°(R d ) subquadratic. With the above regularity result at hand, classical 
semigroup theory |22j allows us to construct a solution to the adjoint problem. 



Proposition 3.6. Let X ^ 0, a G J 

subquadratic. For m > d/2, let tpQ G 
unique mild solution 

cpeC([0,T};L 2 ( 



I with a < 2/(d - 
E"\ V G W m ^(l 



2). 



and U G C°°(R d ) be 
Then, (3.9) admits a 



*))■ 
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Proof. First, we study the homogenous equation d^P(ip(a),a)^ = 0, associated to 
(3.9). It can be written as 

fit* - -iH£ + B(t)£, 

where 

fl(t)£ := -i (X(a + 1)H 2CT C + Aa|V| 2CT - 2 V 2 ? + a(t)V(x)£) . 

The operator —iff : S 2 — > L 2 (M d ) is simply the generator of the Schrodinger group 
S(t) = e- tHt . On the other hand, for any t G [0,T], B(t) is a linear operator 



on the real vector space L 2 (R d ), equipped with the inner product (2.3) (the same 
would not be true if we would consider L 2 (M. d ) as a complex vector space). In 
addition, B(t)* = B{t) is symmetric with respect to this inner product and the 
same is true for iB(t). Since V € W m '°°(R d ), a e L°°(0,T) by assumption and 



ijj e L°°((0,T)xR d ) in view of Lemma 3.4 we infer B e L°°(0,T;£(L 2 (R d ))). The 



operator B(t) may therefore be considered as a (time-dependent) perturbation of 
the generator —iH. 

Following the construction given in Proposition 1.2, Chapter 3 of [22], we obtain 
the existence of a propagator F(t, s), i.e. a family of bounded operators 

{F(t,s) :L 2 (R d )^L 2 (R d )} s . tmT] 

which are strongly continuous in time and satisfy F(t,s) = F(t,r)F(r, s). This 
propagator F(t, s) is implicitly given by 

F(t, s) = e - lH{t - s) + f e~ lH(t - T) B{T)F{T, s) dr 

J S 

and solves the homogeneous linearized equation in the sense that 
^F(t,8)£=(-iH + B(t))F(t,8)t 

weakly in (E 2 )* for every * e L 2 (R d ) and almost every t € [0, T]. Clearly, it 
provides a unique mild solution £(t) = F(t, s)ip(s) of the homogenous equation. 



Duhamel's formula applied to the adjoint problem (3.9) consequently yields 



( , 12) ^.^^^jTw^* 

Under our assumptions on ip and A we have that 



which in view of Duhamel's formula (3.12) implies the existence of a mild solution 
f G C([0, T]; L 2 (M. d )). Uniqueness follows from linearity and the uniqueness of the 
homogeneous equation. □ 

4. Rigorous characterization of critical points 

A classical approach for making the derivation of the adjoint system rigor- 
ous is based on the implicit function theorem. The latter is used to show that 
d^,P{%jj(a), a) is indeed invcrtiblc, but it requires the identification of a linear func- 
tion space X such that 

P : T(0,T) x H^O^T) -> X ;(ip,a) ^ P(ip,a), 

and 

d^P^.a)- 1 :X->T(0,T). 



In other words, we require the solution of (3.9) with a right hand side in X to 
be in T(0,T). It seems, however, that the linearized operator d^P^ja)^ 1 is 



16 M. HINTERMULLER, D. MARAHRENS, P. A. MARKOWICH, AND C. SPARBER 



not sufficiently regularizing to allow for an easy identification of X. Therefore we 
shall not invoke the implicit function theorem but rather calculate the Gateaux- 



derivative J'ia) directly. (We do not prove Frechet-differentiability; see Remark 4.4 



below.) To this end, we shall first show that the solution ip = ip(a) to ( 1.1 1 depends 
Lipschitz-continuously on the control parameter a. This will henceforth be used to 
estimate the error terms appearing in the derivative of J {a). 

4.1. Lipschitz continuity with respect to the control. As a first step towards 
full Lipschitz continuity, we prove local-in-time Lipschitz continuity of ip = ip(a) 
with respect to the control parameter a. 

Proposition 4.1. Let A ^ 0. a £ N with a < 2/(d - 2), and U £ C°°(R d ) 
be subquadratic. For m > d/2, let V € W m ^(R d ) and rp,ip € L°°(0,T;£ m ) be 
two mild solutions to (1.1 1, corresponding to initial data ipo,ipo € £ m and control 
parameters a, a £ H l {Q,T), respectively. Assume that 

l|a|Uj,l|a||irj.ll^(* ) -)l|E'»,IIV'(*.OI|E«. <M 

for some given M ^ 0. Then there exist r = r(M) > and a constant C — 
C(M) < +oo, such that 

(4.1) U - ^|U«(i t;Sm) < C (||#) - +W&- a\\ H }) > 

where I t := [t,t + r] n [0, T]. In particular, the mapping a H ► ip(a) £ T(0,T) is 
continuous with respect to a £ H (0, T). 

Proof. To simplify notation, let us assume t + t ^ T. By construction, there exists 
a t > depending only on M, such that ip\j t is a fixed point of the mapping 

iP^S{-)ip -i J S(-~s)(X\iP( S )\ 2,J iP( S )+a(s)ViP( S )) ds, 

which maps the set 

Y = {^£ L™(I t ; S m ) : |MU~(/ t;S ™) < 2M} 

into itself. Of course, the same holds true for ip and a in place of ip and a, respec- 
tively. In particular, the embedding T, m (R d ) <-> L°°(R d ), m > d/2, yields 

IMIl~(j«xb<) < 2CM. 
To proceed further, we recall the following result, which can be proved along the 
lines of 8, Lemma 4. 10. 2]. 

Lemma 4.2. Let M > 0, a £ N, and m > d/2. Then there exists a constant 
C(M) > 0, such that for all ip, ip £ S m satisfying \\iP\\l™, ||V>IUg° ^ M, it holds 
that 



\iP\ 2 °iP-\iP\ 2 °iP\L m ^C(M)\\iP-iP\\ 



S" 



In other words, ip H ► |^| 2tT ?/> is locally Lipschitz in S m . 

Subtracting the two fixed point expressions for ip and ip gives 
^(s)-V(s) = S(s 

S(s-r)(\(\iP\ 2 ' J iP- \ip\ 2a ip) + V(x)(aip - aip))(r) dr 

for all s £ [t,t + r]. Taking the L°°(I t ; £ m )-norm and recalling Lemma 3.3 together 
with ||V»0OI|e»,||$(*)||e™ < 2M, for s ^ t + r, yields 

\\4> ~ ^IU-</,;E~) < C\\i>{t) - iP(t)\\vm + 2M\\a - a\\ H} \\V\\ wr ,~ 

+ Ct (C{2M) + \\a\\ Hl \\V\\ w ™,~) U - iP\\ L °o it , t+T ., sm) , 
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where C(2M) is the constant appearing in Lemma 4.2 with 2M replacing M. Since 



^ M, the estimate (4.1 1 follows from possibly choosing r even smaller. 
Finally, we show the continuity of the map i? 1 (0,T) — > T(0,T),a i— > ip(a). Set 

f, := inf {0 s$ t < T : limsup ||0 ~ ^IU-(o,t;&») > 0}, 



a— >a 



with the convention inf := +oo. We have to show that t* = +oo. Assuming 
i* < T < +oo, fix M' > M such that ||V'IU«(o 1 T;S"') < M \ let r' = r(M' + 1) > 
be chosen as above, with M' + 1 replacing M . Furthermore let At = r'/2. The 
definition of of t* yields 

limsup ||0 - -0IU~(o,t,-At;S™) = 0. 
In particular, it holds that ||0||L°°(CM,»-At;E m ) ^ M' + l for all (a— a) small enough. 



But now we see that the Lipschitz continuity (4.1 1 is satisfied by tp and ip and such 



controls a, a on the interval [t* — At, t* — At + t'}. Hence 

limsup ||0 - VIU«»(0,t,-At+r';S™) = ° = 

a contradiction to the definition of t*. Hence we must have t* = oo, and continuity 
holds. □ 

As a direct consequence of this continuity result, we obtain uniform boundedness 
of the solution 0(a) on compact sets in a G H X (Q,T). Of course, bounded sets in 
il 1 (0, T) are in general not compact and thus we have to restrict ourselves to finite- 
dimensional subsets. 



Corollary 4.3. Under the assumptions of Proposition 4-1 let 5 a G i? 1 (0,T) with 
<5q(0) = be a direction of change for a and let ip(a + sd a ) be the solution to ( 1.1 ) 
with control a + e5 a and initial data 0o G T, m , m > d/2. Then there exists M < oo 
such that 

\\tp(a + E<5 Q )|| L oo (0!T . Sm) < M, Ve G [-1,1]. 

Remark 4.4. This bound on finite dimensional subsets of i? 1 (0,T) is the reason 
why we can only prove Gateaux-differentiability. If we had a bound on 0(a) in 
the £ m -norm which was uniform in t ^ T and ||a||#i ^5 M, we could prove 
Frechet-diffcrentiability. For our further analysis, however, this will not be of any 
consequence. 

Now we are ready to prove Lipschitz-continuity of the solution 0(a) with respect 
to the control parameter a G i/ 1 (0,T) on the whole control interval [0, T]. 

Proposition 4.5. Let A ^ 0, a G N with a < 2/(d - 2), and U G C°°(R d ) be 
subquadratic. For m > d/2, let V € W m >°°(R d ), O € E"\ a G ff^O.T), and 
■0 = 0(a) G L°°(0,T;i; m ) 6e i/ie solution to (II]). Sei = -0(a) where for any 
e G [—1,1], we let a := a + eS a with S a G i/ 1 (0,T) suc/i that 5 a (0) = 0. ITien i/iere 
exists a constant C > 0, s«c/i £/iaf 

(4.2) ||0 - 0||l~(o,t ; S'") «S C||a - a|| H i (0i T) = C|e|||^llffi(o,T)- 



In other words, the solution to (1.1) depends Lipschitz- continuously on the control 
a for each fixed direction 5 a . 



Proof. Since Corollary 4.3 provides a uniform (in e) bound on ||0||_L°c(o,T;£ m )j the 



quantity r in the local Lipschitz estimate (4.1) is now independent of e and t and 



the estimate indeed holds on every interval [t, t + r], i.e. 

110 - 0llL~(t,t+r;£™) < G (ll^M ~ 0WI|s m + ||5 - a|| H i(t,t+r) 
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Since both solutions ip and ip coincide at t = 0, finite summation of this estimate 
over intervals [nr, (n + 1)t] yields (4.2 1. □ 

4.2. Proof of differentiability and characterization of critical points. We 

are now in a position to state the second main result of this work. 

Theorem 4.6. Let X > 0, a £ N with a < 2/(d - 2), and U £ C°°(R d ) be 
subquadratic. In addition, let tp & £"\ V £ W m -'°°(R d ) for some m £ N } m ^ 2. 
and a 6 Zf^O,! 1 ). 

T/ien i/ie solution of (1.1 1 satisfies ip £ L°°(0, T; £ m ) and i/ie functional Sf(a) is 
Gateaux- differ entiable for all t £ [0, T], wzi/i 

(4.3) J" (a) = Re / Jp(t, x)V{x)^(t, x) dx - 2- (d(t) (72 + 7i^ 2 (*))) , 

J Rd at 

in the sense of distributions, where Oj[t) is the weight factor defined in ( |1.9[ ) and 
ip £ C([0, T]; L 2 (R d )) is the solution of the adjoint equation 

idtf = — ^A^ + U(x)(p + a(t)V(x)ip + X(o~ + l)\ip\ 2cr (p + Xa\i(j\ 2 ' 7 ^ 2 ip 2 Tp 

(4.4) ^ 

+ 4 ll (a(t)) 2 u(t)V(x)4,, 
subject to Cauchy data tp(T,x) = Ai(ip(T, -),Aip(T, -))l 2 Aip(T,x). 
Remark 4.7. When compared to the assumptions of Theorem |2.1| the result of 



Theorem 4.6 requires additional regularity (and stronger decay) of the initial data 
ipo and the potential V (plus, we need to restrict ourselves to a £ N). Note that the 
requirement m £ N and m 2 implies m > d/2 for d = 1,2,3 spatial dimensions. 



Proof. We need to prove that J'ipt) is of the form (4.2 1. For this pur pose , let 

and 



4.5 



ip = ip(ct), ip = ip(ct) with a = a + ed a , satisfy the assumptions of Lemma 
consider the difference of the corresponding objective functionals J {a), J (a). This 
difference can be written as the sum of three terms 

J(a) -J{a) =1 + 11 + 111, 

where we define 

I := $(T), Ai>{T))l 2 - (</>(T), Ai>(T))l 2 , II := 72 f \&(t)) 2 - (a(t)) 2 dt, 

Jo 

and 

III := 71 jf (Mt)) 2 (J V(x)\ip(t,x)\ 2 dxj dt 

-71 / (a(t)) 2 (f V(x)\1>(t,x)\ a dx] dt. 



The general strategy will be to use the Lipschitz property established in Lemma 4.5 
and rewrite the terms I, II, and III in such a way that 

J{6t) — J{a) = linear terms in (a — a) + 0(\\a — ot\\ H i). 

Since a = a + eS a and thus 0(||6j — oj||?ti) = 0(e 2 ), the limit e —> then yields the 
desired functional derivative. 

We start by considering the term I. It can be rewritten in the form 

I = $(T), A^{T)) 2 Ll - (yj(T),A^(T)) 2 Ll 
= 2{^{T),A^{T)) Ll ($(T),A${T))v m - (^(T),^(T)> 



LI 

(^(T), Aip(T)) L 2 - (^(T),Ai;(T)) Ll 
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Using the essential self-adjointness of A, the terms within the parentheses yield 

(i>{T),Ai,{T)) Ll - (i>{T),A^{T)) Ll 

= 2${T) - ip{T), A^{T)) Ll + ${T) - ^{T),A{i>{T) - ^(T)))m. 



Using the Lipschitz-estimate (4.2), we obtain 

$(T) - ^(T), A$(T) - ^{T))) Ll \ < \\A\\ c{11 ,li)U{T) - < C^IIM**, 

and hence 

(i;(T),Ai>(T)) L 2 - {4>(T),A^{T)) Ll = 2(i>(T) - ^(T), A^T))^ + 0(||a- a^i). 

Squaring the above result and plugging it into our expression for / consequently 
yields 

(4.5) l = A{^{T),A^{T)) Ll (i>(T)-i>{T),Ai>{T)) Ll + 0(\\a- af^). 
Next we consider II, which can be written as 

r T r T ■ 2 

II = 2 72 / d(t) (a(t) - a(t)) dt + j 2 / (&(t) - a(t)) dt 



2 72 / a(t)(&(t)-a(t))dt + 0(\\a-a\\ 2 H i). 



The first term in the second line is thereby seen to be of the form given in (4.2). 
Finally we consider III, which in view of definition (1.9) can be written as 

III = 7i f \(&(t)) 2 -(a(t)) 2 )u, 2 (t)dt 
Jo 



7i 



WW 



o 



V{x)\ip{t,x)\ 2 dx -uj 2 (t) dt 



As before, we can expand these terms using quadratic expansions in both ip and 

is 

2 



a. In view of the Lipschitz estimate (4.2), any quadratic error ||-0 — ip\\ 2 LaaL2 
bounded by 0(||<S — a||^i) and hence we obtain 



III 



(4.6) 



4 7l j {a{t)) 2 uj(t) (Re j (($ - ^)V^) (t,x) da?) 



dt 



+ 271 / (&(t) - a(t)) a{t)uj 2 (t) dt + 0(\\ a — all tti ). 



Here the second term on the right hand side is linear in (a — a) and hence of the 
desired form. In order to treat the first term, we note that the expression 

47i (W)) 2 [ V(x)m,x)\ 2 dx] V(xMt,x) 



appears as a source term in the adjoint equation (4.4). Thus we obtain 

rT _ 

)dt 



(4.7) 



4 7l j {a{t)f w(t) (Re J - ^)Vif) (t, x) dx) 

<p(t,x) (d^P(^,a)(^-^)) (t,x) 



Re 

- Re 



JR d 

i7p(T,x)(ijj(T,x) — tp(T,x))dx, 



where we recall that d^P(tp, a) denotes the linearized Schrodinger operator ob- 
tained in (3.5). The last term on the right hand side of (4.7) stems from the 
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boundary condition at t — T. Note that the boundary term at t — vanishes since 
^(0) = Vo = V>(0) by assumption. We recall that <p G C([0, T];L 2 (R d )) and 

€ L°°(0, T; E m ) n W 1,oo (0, T; £ m ~ 2 ), with m > 2, 



and hence the right hand side of (4.7) is well-defined. In addition, since both ip 
and ip solve the nonlinear Schrodinger equation (1.1 1, we can write 

d^P(^, a)$ -if>) = idS ~ i>) ~ - V) ~ V{x){a{t)j) - a(t)i/>) 
(4.8) + A|-0| 2CT -0 - A|^| 2CT ^ + (fi(i) - a(*))V(x)^ + V) 

where the remainder g(tp, tp) is given by 



A 



g(i>, i>) = $\ 2 °i> - IVI 2 ^ - (t + l)l^| 2 "(^ - V) - ^| 2CT - V - 1>) ■ 

Since 2a ^ 2 by assumption, H^Hl^l^ 7 IIV'IIl^ls ^ C m view of Corollary 
and E m =-> L°°(E d ), the remainder can be bounded by 



4.3 



In addition, since <p G C([0, T]; L 2 (M d )) and S m C H m (R d ) 
that 



|¥»(t,a:)||^,a;)-V(t,a;)| 2 da;< |M 



L 4 (M d ), we find 



0(\\a-af H} ). 



Furthermore, the contribution of (a(t) — a(t))V(x)^ in (4.8) equals 

(fi(t) - a(t))V(x)V + (a(t) - a(t))V(a:)(Vi - V), 
where the latter term can be estimated by 0(\\a 
this shows that 



as before. In summary, 



(4.7)=/ (a(t) - a(i))Re / ^(t,a;)V(a;)V>(*, x)da;d* + 0(||a - all^) 
(4.9) — Jo JR d f 

- 4(V(T),^(T)) i j $(T) - V(T), A^j{T)) Ll , 

where we have used the fact that the data of the adjoint problem at t = T is given 
by 

<p(T, x) = U{ip(T, ■), AtP(T, -)) Ll Aip(T, x). 

Thus, we infer that, up to quadratic errors, the second line in (4.9) cancels with the 
terms obtained in (4.5). Collecting all the expressions obtained for I, II, III and 
taking the limit e — > 0, we have shown that S(a) is Gateaux-differentiable with 
derivative J' {a) given by (4.2). This concludes proof of Theorem 4.6 □ 



Equation (4.3) yields the following characterization of the critical points a* G 
H x (0,T), i.e. points where J' (a*) = 0. 



Corollary 4.8. Lettp^, be the solution of (1.1) with control a*. Also, let ip* be the 

corresponding solution of the adjoint equation (4.4), and denote by the function 
defined in (1.9) with ip replaced by ip*- Then a* £ C 2 (0,T) is a classical solution 
of the following ordinary differential equation 



(4.10) 



dt 



1 



Re 



(p*(t,x)V(x)ifj*(t,x) dx, 



subject to a*(0) = cxq, d*(T) = 0. 



Remark 4.9. In the case 71 = this simplifies to the expression used in the physics 
literature; cf. H"4"l. 
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Proof. Le t fi € Cq°(0,T) be a te st function with compact support in (0,T). Then, 

imply that there exists a* <E i? 1 (0,T) such that 



2.1 



and Theorem 



4.6 



Theorem 

i7'(a*) = 0, satisfying (|4.3| in the sense of distributions, i.e 



«•(*)£(<) (7a + Tiw2(*))* 



-Re 
2 ./o 



[i(t)ip*(t, x)V(x)ip*(t, x)dxdt, 



where we have used the fact that the boundary terms at t = and t = T vanish 
due to the compact support of n[t). We shall show that the weak solution a* is in 
fact unique. This can be seen by considering two different al(t),a^(t), satisfying 



aj(0) = a 2 (0) 
solves 



ao- Denoting their difference by /3* = a\ — oq, we have that /3*(t) 



/9,(t)£(t) (72 + liu*{t)) dt = 0, for all M G C °°(0,T). 

Since 72 > and 71 0, this implies that /?* (t) — in the sense of distributions. 
However, since a\,al G H 1 ^^) C(0,T), we conclude that & G C(0,T) and 
thus /3*(t) = const for all t G [0, T]. Since /3*(0) = by assumption, we infer unique- 
ness of the weak solution a*(t). On the other hand, standard arguments imply that 



(4.10 ) admits a unique classical solution a* G C 2 (0, T), provided G C 1 (0, T) and 
the (source term on the) right hand side is continuous in time. The latter is ob- 
viously true in view of Proposition |2.3| and Proposition |3.6| In addition, since 
V G W^°°(R d ), we infer that for all ip(t) G E it holds that X (t) ■= (V(x)ip(t)) G E. 



From Proposition 2.3 it follows that 



w,(t) = 2Re / V(x)d t ip*(t,x)ip*(t,x)dx = 2(x(t),'>p*(t)) 



< 



Thus, uj(t) G C 1 (0,T), yielding the existence of a unique classical solution a* G 
C 2 (0, T). We therefore conclude that the unique weak solution a* obtained above 
is in fact a classical solution, satisfying (4.10) subject to a*(0) = cto, d*(T) = 0. □ 



We call a* G H l (Q. T) a critical or stationary point of the problem 
(4.11) min J(a) over a € ff 1 (0, T) 



if J' (u*) = 0, where J7 7 is given in Theorem 4.6 In order the check computationally 
whether a* is critical, one needs to solve (JTTTJ for a — a* to obtain ip* and then 
the adjoint equation ( |4.4| with ip = and a = a* to compute y». Inserting 
(a,ip,(p) — (a*, f/'*, V 3 *) m (4.3) yields J7'(a*) which has to vanish for a* to be 
critical, i.e., (14.101) is satisfied. We therefore call dl.lh , (|4.4h and J4.101 the first 



order optimality conditions associated with (4.11) 



5. Numerical simulation of the optimal control problem 

For our numerical treatment we simplify to the case d = o~ = 1. In this case, the 
first order optimality conditions for our optimal control problem are given by: 



dt 



(d(t) (72+7i^ 2 (<))) 



1 



Re 



ip(t, x)V(x)ip(t, x) dx, 



id t ip + -<9 2 ^ = (17(a) + a(t)V(x))ip + A|^|V, 

id t cp + -d%tp = (U(x) + a(t)V(x))ip + 2X\ip\ 2 if + \4> 2 Tp + 4 7l (d) 2 oj(t)V(x)4>, 

subject to the following conditions: a(0) = cto, d(T) = 0, and 

iP(0,x) =Mx), f(T,x) = ii(iP(T,-),A^(T,-)) L2 M(T,x). 
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In our numerical simulations, the resulting Cauchy problems for Schrodinger- 
type equations are solved by a time- splitting spectral method of second order (Strang- 
splitting), as can be found in [Tj. This computational approach is unconditionally 
stable and allows for spectral accuracy in the resolution of the wave function ^(i, x). 
This is needed due to the highly oscillatory nature of solutions to (nonlinear) 
Schrodinger-type equations. We consequently perform our simulations on a nu- 
merical domain OcK, equipped with periodic boundary conditions. The trapping 
potential U (x) is thereby chosen such that the "effective" (i.e. the numerically rele- 
vant) support of the wave function tp(t,x) stays away from the boundary. In doing 
so, the boundary conditions do not significantly influence our results. A good test 
of the accuracy of our numerical code is given by the fact that the Gross-Pit aevskii 
equation conserves the physical mass (i.e. the L 2 -norm of ip(t)). Indeed, in all our 




numerical examples presented in Section 5.3 below, we find that the L 2 -norm is 
numerically preserved up to relative errors of the order 10~ 13 . 

5.1. Gradient-related descent method. Once a suitable solver for the state and 
the adjoint equations is at hand, our gradient-related descent scheme operates as 
follows. We determine a sequence of descent directions C i? 1 (0,T), i.e., for 
every k G N 

J(a k + 8 k a ) < J(a k ) ee J{ip{a k ), a k ). 

is satisfied. Note that a simple Taylor expansion of J around a k shows that 
(i7'(afc), (^) < is sufficient for 5* to be a descent direction for J at a k . We 
are in particular interested in gradient-related descent directions which satisfy 

MS* = -J'{a k ) 

with M a suitably chosen positive definite operator. 

A rather straightforward choice of M is given by M = d^J(tp, a k ). In this case 
8% is obtained as the solution of the following ordinary differential equation (of 
second order): 



J'(a k ) =2j t (sl(t) f 72+71 Uv(x)\i> k (t,x)\ 2 dx 



with 6^(0) = and <J*(T) = 0. Here ip k (t,x) denotes the solution of the Gross- 
Pitaevskii equation with a(t) = a k {t). With this choice of a descent direction, we 
then perform a line search in order to decide on the length of the step taken along 
6*- In fact, we seek for v k > such that 

(5.1) J{a k + v k 5*) < J(a k )+»v k {j'(a k ),5*) 

with some fixed fj, € (0,1). Within each line search, we determine v k iteratively 
by a backtracking strategy. Thus, the whole procedure amounts to an Armijo line 
search method with backtracking. Of course, more elaborate strategies based on 
interpolation or alternative line search criteria are possible; see, e.g., |23] for more 
details. 

We stop the gradient descent method whenever 

(5.2) \\J'{a k )\\ H -x <TO-L-\\J'{a x )\\ H ^ 

is satisfied for the first time. Here, TOL € (0, 1) is a given stopping tolerance and 
ct\ G i/ 1 (0,T) is the initial guess satisfying the boundary conditions ai(0) = qq 
and ai(T) — 0. As a safeguard, also an upper bound on the number of iterations 
is implemented. 

In our tests, we observe the usual behavior of steepest descent type algorithms, 
i.e., the method exhibits rather fast progress towards a stationary point in early 
iterations, but then suffers from scaling effects reducing the convergence speed. 
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Therefore, often the maximum number of iterations is reached. Thus, we connect 
the first-order, gradient method to a Newton-type method which relies on second 
derivatives or approximations thereof. 

5.2. Newton method. The majority of iterations within our simulations are per- 
formed via a second order method, Newton's method, for which we use the full 
Hessian 

M := D 2 a J{a k ) : H\0,T) x fl^O.T) -> K, 
or a sufficiently close positive definite approximation thereof. Note that we can 
also consider the Hessian as a map D 2 a J : if 1 (0,T) — s- ^{Q.T)* . Recall that the 
gradient-related method above simply uses M = d 2 J(ip 7 a k ). 

We derive P> 2 a J formally form the Lagrangian formulation; see Remark 3.1 The 
Lagrangian is given by 

L(tp,a,(p) = J{ip,a) - (<p,P(ip,a)) L 2 , 



where ip is the solution to the adjoint equation (4.4) and P(^, a) is the Gross- 



Pitaevksii operator written in abstract form. Proceeding formally, we find 

hi. 



((DlJ)8 a ,S a ) L 2 = ((d 2 L)5ji,,5j(,} L 2 + ((^ a I)(J tt ,^) 



+ ((d a ^L)6 a ,S^) L 2 + ((d^L)S a ,S a ) L 2, 

where 8^ and 5^ solve the linearized Gross-Pitaevksii equation with controls 8 a , S a 
respectively. In view of the derivation given in Section |3.1| we have 

8^ = i/>'(a)8 a = —d^P(tp(a) , a)' 1 d a P(^(a) , a)6 a , 

and analogously for 5^. Hence we conclude that 

((D 2 a J)S a ,~5 a ) L 2 = ((dlJ)il>'(a)5 a ,i)'(a)6 a ) L 2 + ({d a ^J)8 a ^'(a)8 a ) L 2 

/ rr o\ t .x t .x 

+ ((d^ a J)^(a)6 a ,S a ) L 2 + ((d 2 a J)8 a ,8 a ) L 2 - ( v , (D 2 a P(^j, a)6 a )6 a ) L , ^ 

where 

(D 2 a P(^,a)S a )S a = (8%Py>,a)y/(a)8 a ))(l/(a)8 a ) 

+ (d a ^P(i/j,a)5 a ){tp'(a)S a ) + (d tf)a P('ip,a)(i{j'(a)S Q ))S a , 



since d^,P(tp, a) — 0. All of the terms appearing on the right hand side of (5.3) can 
be evaluated by replacing ip'(a)S a by —d^,P(ip,a)~ 1 d a P(tp,a)5 a . Consequently 
for calculating the action of the Hessian this requires to solve several linearized 
Schrodinger-type equations with different source terms and boundary data. For 
example, the term involving {d a ipJ) can be evaluated by using 

X ■= d. 4 ,P{i>,ay*((d a4 ,J)5 a ), 

which solves the following Cauchy problem 

idtX +\dlx = U (x) X + a{t)V(x) X + 2A|Vfx + X^ 2 X + &nh(t, x)i>, 

where h(t.x) := uj(t)a(t)S a (t)V(x) and 

6 2 J(^,a) 
X{ ' ' 8^(T,x)8a(T) 

Bearing this in mind, we have to solve the following equation for 8* € _ff 1 (0,T): 

(5.4) M8 k a = D 2 a JS k a = -J'(a k ) E H\Q,T)*. 

Hence, we need to invert D 2 J, which, in view of ( |5.3| ) is not directly possible. 
Rather we resort to an iterative method, the preconditioned MINRES algorithm, 
see [H], with the preconditioner {d 2 a J , a))- 1 : H 1 (0,T)* -> H 1 ^^). 
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We emphasize that here we aim to study the behavior of solutions of our control 
problem rather than at optimizing the respective solution algorithm or its imple- 
mentation. 

5.3. Numerical examples. In all our examples, we choose the numerical domain 
Q = [-L, L] with L = 20 and periodic boundary conditions. The number of spatial 
grid points is N = 256. In addition, we set the final control time to be T = 10, 
and we use M — 1024 equidistant time steps. In order to avoid the influence of 
the boundary, we choose a trapping potential U{x) — 30 (^) . The initial guess 
for the control is taken to be just a.\ = in the linear case (A = 0), whereas 
each algorithm in the nonlinear case (A ^ 0) is started from the control obtained 
by solving the linear problem. In our tests of the first-order gradient method, we 



choose TOL = 10 in the terminating condition (5.2) for the whole algorithm, 
/i = 10~ 3 , and a maximum number of 20000 iterations. For the Newton method, 
we likewise set TOL = 10~ 8 and we stop the algorithm after at most 45 Newton 
steps. 

5.3.1. Example: shifting a linear wave packet. For validation purposes, we consider 
the time-evolution of a linear wave packet, i.e. A = 0, whose center of mass we aim 
to shift towards a prescribed point yi G [—L,L]. For this purpose consider a control 
potential 

^ro + lx^ < VseH.L], 

and the observable 

In this case, we find that the algorithm converges well even if we only invoke the 
first order gradient method. Indeed, as we decrease the regularization parameters 
7i,72 <C 1, we approach an optimal solution which, as it seems, cannot be im- 
proved upon. This optimal solution, or, more precisely, its spatial density p = , 
is depicted in Figure [l] (right plot) , where we denote by "target" the function pro- 
portional to 1 — A(x) with k = 0.07 and y\ — —2L/8, such that it has the same 
L 2 -norm as i/V The left plot shows the associated control. 



Figure 1 . Shifting a linear wave paket 



Since this solution seems optimal, the choice of 71 , 72 becomes negligible below 
a certain threshold. Thus, it suffices to consider 71 = and only include the cost 
term proportional to 72. Similar results hold for any other given point y\ € f2, 
provided y\ stays sufficiently far away from the boundary. 
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5.3.2. Example: splitting a linear wave paket. We still consider the linear case, i.e., 
A = 0, and aim to split a given initial wave packet into two separate packets centered 
around yi and y 2 ? respectively. The control potential is chosen as 

V{x) = e - 8 * 2 /i 2 > 0, 

and the observable 



In the following we fix k = 0.07, y\ = —2L/8, and y% — 2L/8. In this case we 
find that the residual of the first order gradient method does not drop below the 



tolerance given in (5.2) before the maximum number of iterations is reached. With 
the Newton method, however, we find a (local) minimum of the objective functional 
J{ijj, a) in less than 20 Newton iterations. Of course there is no guarantee that this 
is a global minimum. 

In order to illustrate our results, we consider the case where 71 = 0, 72 = 
1.5 x 10~ 6 . At the final control time T = 10 we then obtain: 

(A^(T),tJj{T)) 2 L 2 « 2.261 x 10" 3 . 

The spatial density p = \ip\ 2 of the corresponding solution is shown in the right 
plot of Figure [2] The associated control is depicted in the left plot. If, instead, we 



Figure 2. Splitting a linear wave paket with 71 = 



choose 71 = 4x 10 



72 = 1 x 10 9 , we find 
(Aip(T),i/>(T)} 2 L 2 ~ 2.269 x 10" 



and the corresponding solution is given in Figure [3| Here the intermediate state is 
a plot of p{t) at t = 4 = 0.4 x T. 

A direct comparison of the (spatial densities of the) resulting wave functions and 
the respective controls is given in Figure |4j We see that the spatial densities are 
nearly identical, but the variability of the respective control parameters is not the 
same. This is, of course, related to time-evolution of the weight factor oj(t), defined 
in ( 1.9 1, which is shown in Figure [5] for the case of 71 = 4 x 10 -5 and 72 = 1 x 10~ 9 . 
By construction, the time-integral of ui(t) can be interpreted as the physical work 
performed during the control process. We find that compared to the case 71 > 0, 
the term ||-B(-)|||2 is around 30% larger (64.5 versus 49.1) and ||J5(')||| a is around 
twice as large (95.0 versus 43.4) in the case where 71 = 0, yielding a significant 
advantage of our control cost over terms considering the 7J 1 -norm only; see [TJ] for 
the latter. 
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Figure 3. Splitting a linear wave paket with 71 > 




Figure 4. Direct comparison between results 




FIGURE 5. The weight factor uj = J V\ip\ 2 dx over time 

Finally, Figure [6] shows an example of the evolution of the objective functional 
J{ijj, a) over the number of iterations of the Newton method, here for the case 
where 71 = 0. 

5.3.3. Example: splitting a Bose-Einstein condensate. We consider the same situ- 
ation as in the previous example, but with an additional (cubic) nonlincarity. More 
precisely, we choose A = 8 > 0. It turns out that the conclusions are similar to 
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Figure 6. Value of J(V>, a) over number of iterations 

the ones found in the linear case (A = 0). Qualitatively, the main difference is that 
during the time-evolution, the wave function spreads out more because of the addi- 
tionally repulsive (defocusing) nonlinearity. In the linear case, the widest extension 
of the wave packet is always comparable to its final value. Choosing as before 
71 = 4 x 10~ 5 and 72 = 1 x 10 -9 , we obtain the solution depicted in the right plot 
of Figure [7J where we show the spatial density at the times t = 0, t = T = 10 and 
at the intermediate time t = 4. The control is shown in the left plot. In comparison 
to the linear case (A = 0), the observable term in the objective functional J(ip, a) 
is found to be slightly larger. Indeed, we obtain 

(A^(T),ip(T)) 2 L 2 « 3.720 x 10~ 3 . 

This seems to indicate that nonlinear effects counteract the influence of the control 
potential. 




Figure 7. Splitting a condensate with 71 > 

We again compare the present case with the one where 71 = (i.e. no cost term 
proportional to the physical work) and 72 = 1.5 x 10~ 6 . First, we find that 

{Ail>(T),ip(T)) 2 r2 « 3.382 x 10~ 3 . 

a; 

Moreover, ||-E||^2 is about 150% larger (172.1 versus 68.5) than in the case where 
71 7^ 0. Similarly, the total energy ||-E||^2 is around 15% larger (91.8 versus 79.5). 
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5.3.4. Example: splitting an attractive Bose-Einstein condensate. Our numerical 
method allows us to go beyond the rigorous mathematical theory developed in 
the early chapters. In particular we may try to control the behavior of attractive 
condensates, which are modeled by (1.1) with A < 0, i.e. a focusing nonlinearity. 
Here we choose A = — 1, whereas the parameters 71 = 4 x 10~ 5 , 72 = 1 x 10 -9 
are the same as before. The results are shown in Figure [8] (control in the left plot 
and the state at times t = 0, 10,4 in the right plot). The observable part of the 
objective functional satisfies 

(AV>CT),V>(T))| a « 2.143 x 1(T 3 . 



Figure 8. Splitting a focusing condensate with 71 > 



In comparison to the case of a repulsive (defocusing) nonlinearity the final value 
for the observable term (Atp(T),ip(T))^ 2 is much smaller, confirming the basic 
intuition that an attractive condensate does not tend to spread out as much as in 
the repulsive case. 

6. Concluding remarks 

In this work, we have introduced a rigorous mathematical framework for optimal 
quantum control of linear and nonlinear Schrodinger equations of Gross-Pitaevskii 
type. We remark that in the physics literature, L 2 (M. d ) is usually considered as a 
complex Hilbert space, equipped with the inner product (<p,£) = Jm d f{x)^{x)dx, 
whereas we consider L 2 {W 1 ) as a real Hilbert space (of complex functions), equipped 



with (2.3). Note, however, that the expectation value of any physical observable A 



and thus also J(ip, a) is the same for both choices. 

Let us briefly discuss possible generalizations for which our results remain valid. 
First, we point out that in our analysis above, we did not take advantage of the fact 
that 71 > and hence all of our results remain true in the case 71 = 0. However, 
Example 1 5 . 3 . 2 1 shows a significant quantitative difference in the behavior of the cost 
functionals with and without the term proportional to 71. 

Second, it is straightforward to extend our analysis to the case of several control 
parameters, i.e. 

K 

V(t,x) = Y,<*k(t)Vk(x), K>2. 

k=l 

Clearly, for V k € W m "°°{W i ), m ^ 2 > d/2, all of our results remain valid. In 
addition, it is not difficult to extend our framework to cases of more general control 
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potentials V(a(t), x), not necessarily given in the form of a product. Such potentials 
are of physical significance; see cf. [14]. From the mathematical point of view, all 
of our results still apply provided that 

Wa, -)|| vie- < Ci, \\d s a V(a, .)|U- < C 2 , V|s| < 2. 

Note that in this case, the cost term in J(ip, a), which is proportional to the physical 
work performed throughout the control process, reads 

f (E{t)) 2 dt= f {a(t)) 2 (f d a V(a(t),x)\^(t,x)\ 2 dx] dt. 

Jo JO \JR d J 

It is more problematic to provide a rigorous mathematical framework for control 
potentials V(a,x) which are unbounded with respect to x £ Mr. Only in the case 
where V(a,x) is subquadratic with respect to x and in L°°(M. d ) with respect to a, 
existence of a minimizer can be proved along the lines of the proof of Theorem |2.1| 
More general unbounded control potentials V(a,x) definitely require new mathe- 
matical techniques. Note that in this case, even the existence of solutions to the 
nonlinear Schrodinger equation is not obvious. 

Finally, we want to mention that it is possible to extend our results (with some 
technical effort) to the case of focusing nonlinearities, A < 0, provided a < 2/d. 
The latter prohibits the appearance of finite-time blow-up in the dynamics of the 
Gross-Pitaevskii equation. Clearly, the optimal control problem ceases to make 
sense if the solution to the underlying partial differential equation no longer exists. 
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